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We propose a new method of studying a real-time canonical evolution of field- 
theoretic systems with boundary coupling to a realistic heat bath. In the free-field 
case the method is equivalent to an infinite extension of the system beyond the 
boundary, while in the interacting case the extension of the system is done in linear 
approximation. We use this technique to study kink-antikink dynamics in 99^ field 
theory in 1-|-1 dimensions. 
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I. INTRODUCTION 



There is a growing interest in transitions over energy barriers in field theories at fi- 
nite temperature. The motivation comes from high-energy physics and cosmology {e.g., 
baryon-number violating sphaleron transitions in electroweak theory 0), as well as 
from condensed-matter physics {e.g., current-reducing fluctuations in one-dimensional su- 
perconductors [0). Such thermally induced transitions usually involve collective excitations 
like kink-antikink pairs or sphalerons, and therefore occur in a nonlinear setting. The non- 
linearity, in turn, often renders analytic techniques useless, forcing one to resort to numerical 
(lattice) methods. On the other hand, many questions of interest can be answered by con- 
sidering the classical regime in which the temperature is much larger than the energies of 
the field quanta at the scale of a relevant collective excitation [^. 

A fundamental quantity of interest is the corresponding transition rate. Unlike static 
observables, it cannot be determined from a given canonical ensemble of field configurations. 
Instead, one must follow a real-time evolution of the field-theoretic system in the heat bath. 
Obviously, the rate would depend on the properties of the heat bath and its coupling to the 
system. This was indeed confirmed by recent numerical studies in which the heat 
bath was implemented through Langevin equation. In most naturally occurring situations 
a system interacts with its environment through the boundaries, and the environment is an 
infinite extension of the system itself. This consideration dictates our choice of a heat bath 
in this work. Ideally, the heat bath we model should possess two main properties: (a) waves 
traveling through the boundary out of the system should be completely absorbed, as if there 
were no boundary and the system had an infinite extension; and ( b) the waves traveling into 
the system should be thermally distributed. Our construction of a heat bath ensures these 
two properties in the free-field case, while for an interacting field they hold in the linear 
approximation. The similarity between our linearized heat bath and the true extension of 
the interacting system becomes better with decreasing temperature. Note that the standard 
Langevin dynamics, imposed either in the bulk of the system or at the boundary [|ri|] , will 
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not lead to the properties a and b. 

While the idea of mimicking a natural heat bath is very generally applicable, we restrict 
ourselves in this article to the case of a scalar field in one spatial dimension. In Section II we 
derive the boundary force exerted by the heat bath on this one- dimensional system. Next, 
we show in Section III how our construction of the realistic heat bath can be implemented 
numerically. As an example, in Section IV we apply our technique to study kink-antikink pair 
nucleation in theory, a subject that has been extensively investigated by other methods 
10]- 1^. We summarize and discuss our results in Section V. 



II. THE METHOD 

Consider a lattice system consisting of a (possibly self- interacting) scalar field (/?„, with 
a minimum of its potential V{(p) ai if = v and mass m corresponding to that minimum. For 
definiteness, let it reside on the n < sites of a one-dimensional chain, and be coupled to a 
heat bath at n = 0. Let the heat bath be a free massive scalar field system on the positive 
n sites with equation of motion (the lattice spacing is a) 

5 h m = (1) 

for n > 1 and boundary condition 

Mt) = Mt)-v = fit). (2) 

Given the solution $n(^), there will be a reaction force exerted by the heat bath on the 
system: 

F{t) = ($1 - $o)/al (3) 

The $0 contribution to F{t) is a harmonic force. The $i contribution, called /mem (the 
memory force) in the following, is, as we shall see immediately, a response of a more general 
nature. As the heat bath is a free field system, this response should be linear. Furthermore, 
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it should causally depend on f(t). In other words, there exists a response function 
{t > 0) such that 

/mem(t)= ^ f it')xit - t')dt' . (4) 

J —CO 

In order to determine xit) l^t us take f{t) = e"*. It follows from Eq. |I| that $n(^) = 
exp(st — kgan), with 

kg = — arcsinh f ^a/s^ + . (5) 



a V2 

Using Eq. |^ we then find that should satisfy 



1 /■* 

a' 



^6^'-'^'^= e'''xit-t')dt' = e'\Cx){s), (6) 



where £ denotes the Laplace transform. It follows that 



^C-' (i^l + y-^y) with y = a'is' + m')/A. (7) 



The integral in Eq. ^ is most easily computed in Fourier space. We define ft{uj) = 
/-oo *''/(^')'^^' xi^) = I^oo^^'^^xi't)dt^ where we are free to take any extension 
of x{t) for t < 0. Then it follows from the convolution theorem that 

/mem(t) = ^ f X{^) ft{^^)duJ . (8) 
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For numerical purposes it is best to extend x(t) to t < such that x(^) vanishes outside 
a finite interval. To this end we choose = — x(— t), so x(^) becomes purely imaginary 
(essentially the sine transform). The latter can be obtained immediately from the Laplace 
transform (Eq. |^) by the 7r/2 rotation of s in the complex plane; taking the imaginary part 
then gives 

x{uj) = isign{uj)\J (tu^ — m'^){l + a^{m? — c<J^)/4) (9) 



for frequencies m < |ct;| < \Jiv? ^ ^ija^ corresponding to propagating modes and x(ll') = 
outside this range. 



We expect the heat bath response to be finite for any bounded f(t) (Eq. ^). If so, x(^) 
must decay sufficiently rapidly in the distant future. This indeed is the case: using Eq. |^ 
and performing the inverse Fourier transform we find in the saddle-point approximation 
x{t) oc sin(mt + tt/4) as t ^ oo. 

The dissipation by the heat bath can equivalently be described by a linear response to 
the momentum ttq at the boundary, rather than the field. Explicitly, it follows from Eq. ^ 
that 

/mem(t) = " ^{t - t')Tl^{t')dt\ (lO) 

J —CO 

where il){t) = J^xi'^^ldt' or, equivalently, {Ctjj){s) = s~^{Cx)is), so that 

i^{uj) = tuj-'x{oo). (11) 

We now turn to the second contribution to the boundary force, the random force /ran- 
It is a Gaussian random variable whose properties are defined by its time autocorrelation 
= {fran(t)fran(t + 1")). Thcsc are most easily determined through the fiuctuation- 
dissipation theorem, stating 

Cir) = d^{\T\) (12) 

{6 denotes the temperature). Numerical implementation of /ran amounts to generating 
Gaussian noise with this autocorrelation function. 

One can check explicitly that the time correlation function of the total boundary force 
F + /ran IS cqual to the average of (($„(0) — $„_i(0))(($„(t) — $„,_i(t))/a'^ for a canonical 
ensemble, calculated for an infinitely extended free field. This is to be expected, as ($„ — 
^n-i)/0''^ represents the force between neighboring sites. 

We conclude this section by discussing the continuum limit of the heat bath response. 
Obviously, F(t) diverges as a — > 0. Instead, aF(t) is a well-behaved quantity whose limit 
is simply dx^{x = 0,t), where x is the continuum spatial coordinate. In other words, if we 
prescribe the boundary field motion, the heat bath response will determine the boundary 
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field spatial derivative in such a way that the waves traveling out of the system will not 
be reflected at the boundary. The corresponding response function is found following steps 
analogous to Eq. ^Eq. ^ We write down the continuum equations of motion 

$ - + 771^$ = > 0) (13) 

together with the boundary condition $(0,t) = f{t) and require 

a.$(o,t)= /* fit')xcit-t')dt'. (14) 

Taking f{t) = exp(st) then yields Xcit) = C^^{\/ s"^ + m?), or Xcif) = 5' {t—Q^)+mt^^ Ji{mt) . 
It is an easy exercise to verify that the same result is obtained by first computing the response 
function for the lattice boundary field derivative and then taking the limit a — > 0. 

III. NUMERICAL IMPLEMENTATION 

Our numerical implementation of the boundary force is as follows. The memory force 
/mem(^) is computed using Eq. ^, that is, in Fourier space. The function x(^) is given by 
Eq. ^. In practice, the integral Eq. ^ is replaced by a finite sum over discrete values of uj, 
separated by an increment Acj. As a result, the response function x(t) becomes periodic 
in t with a period lnjlS.LO. Beyond t = 271/ Au our approximation of x(t) is incorrect, and 
the lower limit of integration in Eq. ^ should be cut off. We achieve this by changing the 
definition of ft{(^) to 

kco) = f e-(*'-*)/(t')rft', (15) 

where T < l-n j Alo. At the same time, T should be large enough so the discarded part of the 
integral Eq. ^ is negligible. With the new definition, /((cj) satisfies the equation of motion 

Ui^) = -tuftiuj) + fit) - e-'-^fit - T). (16) 

The random force /ran is computed by convolving white noise of unit power spectrum with 
a function R{t) whose Fourier image is given by \Jc{ijj), where C is the Fourier transform 
of C{t) given by Eq. |12| [Q. This can be done very efficiently using FFT algorithms. 
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With both ingredients of the boundary force in place, we can now write down and 
integrate the equations of motion for the one- dimensional system immersed in the heat 
bath. These equations have the standard form 

V^n 5 V (V9„,) = (17) 

in the bulk of the system, but should be modified at the boundaries. For example, at the 
left boundary we have 

V^O — — — ^'(V^o) = /mem,left + /ran,left, (18) 

and the right-boundary analog is obvious. The system of equations of motion for the field is 
supplemented by Eq. ITB governing the evolution of the memory forces at both boundaries. 



We integrate equations |T^, |T8|, and |T6| using the second-order Runge-Kutta algorithm. 
While this way of updating ft{uj) is efficient compared to computation of the full integral 
Eq. it is not accurate enough to maintain correct phases of ft{^^) for times much longer 
than 27t/Auj. To ensure stability, we therefore adjust the values of ft{^^) by computing 
Eq. |1^ every T time units. 

We have tested the action of /mem by evolving an initially hot system with /ran omitted. 
This corresponds to cooling in the zero-temperature heat bath. The resulting evolution is 
very similar to that of a benchmark run where we substitute for the heat bath a real cold 
free-field system with large volume. This similarity holds separately for every momentum 
mode. The method allows cooling the system to at least 10"^ times its original temperature. 
Cooling curves for a system with a real and with a simulated heat bath are shown in Figure 

B 

We then turned /ran on and tested it by heating a cold field configuration to a prescribed 
temperature. Again it compared well, mode by mode, with a benchmark run, in which we 
coupled a cold system to a large real heat bath whose temperature the system eventually 
reached (let us stress that using the real heat bath is much more costly in terms of CPU 
time). Both heating curves are shown in Figure ^. 



Finally, we have shown numerically the self-consistency of our method by comparing the 



motion of the endpoint field to that in the middle of the system . As Figure ^ shows, 
the autocorrelation curves of the two fields are very close to each other, meaning that the 
simulated free-field heat bath closely approximates the real one. 



IV. KINK-ANTIKINK DYNAMICS IN 1+1 DIMENSIONS 

We have applied our method to investigate kink-antikink dynamics of (p^ theory whose 
Lagrangian is 

L = 1(^^-1)^) (19) 

in suitably chosen units 0. Kinks and antikinks are finite-energy solutions of the equations 
of motion interpolating between the vacuum values of the field ip = ±1. Explicit func- 
tional form of the static kink is v^±(x) = tanh(±^). The spatial extension of a kink is 
approximately and the kink mass is M = ^8/9. 

Following 1^, we chose a system of = 400 sites with lattice spacing a = 0.5. The 
simulations were performed at seven values of inverse temperature (3 between 3.0 and 6.0. 
Our Runge-Kutta time step was 0.005. For each value of the temperature we started from 
an ordered system and warmed it up for 2500 time units followed by 2 x 10^ time units over 
which we measured the number of kinks n. For the latter we used the same definition as in 
the number of zeros in a field configuration smoothened over the physical distance of 
AL = 5. 

The average kink number (n) at a given temperature should not depend on the properties 
of a heat bath. Our measurements of this observable, shown in Figure ^ are indeed close 
to those of 0]. Note, however, that our measurement errors are much larger, especially for 
low temperatures, even though our data sample is as big as in [^]. We estimated the errors 
using a jackknife method with the block size varying over a very long range; we are therefore 
confident that the autocorrelation of our data is properly taken into account. Moreover, we 
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studied a microcanonical evolution of our system at the energy roughly corresponding to 
(3 = 4.5, with the error estimate similar to that of the corresponding canonical case. 

The temperature dependence of (n) may be interpreted in terms of the effective kink 
mass. Namely, one expects M 



{n) oc ^/3exp(-/5Meff). (20) 

It was found in earlier work that Mcs < M. It was also indicated that M^s is temperature 
dependent 0. Both features find further evidence in our study. If we try to fit all our 
measurements of (n) to Eq. |20| at once, an unacceptably low goodness-of-fit results. The 
situation improves dramatically if we exclude the highest-temperature point from the fit. 
We then find M^s = 0.695 ± 0.0095, or M^s = (0.737 ± 0.010)M, in good agreement with 
Alternatively, we can use pairs of consecutive values of (n) to extract Mcs- The result, 
presented in Figure |], shows the tendency of Mes to decrease at higher temperatures, in 
agreement with findings of Q], [^]. 

Another interesting quantity we extract from the kink-antikink number time history is 
a kink lifetime, i.e. the autocorrelation time r of n. The latter is usually obtained by fitting 
the n autocorrelation function {{n(t) — {n)){n{0) — (n))) to a single exponential of the form 
exp(— t/r). If the time history exhibits more than one time scale (as is the case for n), r 
can only be given an average, or effective meaning. The existence of such multiple scales 
also makes a single-exponential fit to the autocorrelation function extremely difficult. A 
multiexponential fit to noisy data is not a practical possibility. We use an alternative way of 
determining r, closely related to the integral definition of the autocorrelation time. Namely, 
if At is a time interval between two consecutive measurements of n, we expect 

N-l 

w{N) = {{Y. - N{n)f) (21) 

^{{n-{n)f) ^ 



_tanh(f) 2sinh^(f) 
Obviously, for large w{N) approaches a random- walk behavior. For a given A^ we de- 
termine w{N) from our data set and solve Eq. ^ for r. For t = NAt » r the result 



is approximately independent of A^, and we take it as an estimate of the kink lifetime. A 
typical dependence of r on t is shown in Figure |^. Note that large values of t for which the 
plateau is reached indicate the existence of multiple time scales in the kink-antikink number 
fluctuations. 

Unlike the exponential fit, this method also allows a well-defined error-estimating proce- 
dure for T. In particular, we apply the jackknife technique. Note that the lag values t used 
to determine r are of the order of 5000. This is to be compared to our n time history length 
of 2 X 10^. We therefore only have a small effective number of independent measurements of 
r, and our error estimate cannot be very accurate. Conservatively we can expect the errors 
of T to be correct within a factor of 2. This might explain their inhomogeneous dependence 
on the temperature. 

Kink-antikink pair nucleation can be viewed as a multidimensional analog of a particle 
escape over a barrier 0, [§], with the kink lifetime related to to the effective barrier height 
B i: 

T oc exp {(3U) , (22) 

where U = B - M^s. From our data (Figure ^ we find U = 1.01 ± 0.06 = (1.07 ± 0.06)M, 
slightly higher than U = (0.85±0.15)M of [Q, obtained by solving a low-viscosity Langevin 
equation. Both that value and ours are inconsistent with the naively expected B = 2M. 
More work is required to explain this discrepancy ||12|. More importantly, however, our 
results show no exponential suppression of the kink-antikink pairs nucleation rate with 
growing temperature, in agreement with analytical predictions and earlier numerical work 



V. CONCLUSIONS AND OUTLOOK 

In this article we presented a method of modeling naturally occurring heat baths. While 
our presentation concentrated on a scalar field in one spatial dimension, the principles un- 
derlying our construction of a heat bath do not depend on the dimensionality or the field 
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content of a theory in question. In any case, one can determine the memory force exerted by 
a heat bath by studying the hnear response of the latter to the field motion at the system 
boundary. The corresponding random component of the force may then be found using 
fluctuation-dissipation theorem. The only new feature to appear in dimensions higher than 
one is related to the connectedness of the boundary: the field motion at different points 
of the boundary will be correlated in a way consistent with causality. This is, however, 
a technical difficulty, not a conceptual one. The work on extending our method to other 
systems is currently in progress. 

We have verified numerically that our simulated heat bath thermalizes correctly both 
linear and nonlinear systems. As an application, we considered the dynamics of kink-antikink 
pairs in (/?^ theory. Our measurements of the kink density agree well with those obtained 
by solving Langevin equation, as one would expect for an equilibrium quantity independent 
of a heat bath implementation. The kink lifetimes we measure are close to those following 
from the low-viscosity Langevin dynamics. This is again to be expected for sufficiently large 
systems: as the system size grows, the influence of the boundary heat bath on the dynamics 
decreases. The Langevin analog would then be decreasing viscosity. 

To conclude, we emphasize again an important advantage of our simulated heat bath 
over the Langevin method (including its zero- viscosity microcanonical limit): our heat bath 
is not arbitrarily chosen and involves no free parameters like viscosity. Rather, it is, to the 
best of our knowledge, the closest known approximation to a natural situation, in which 
open systems are immersed in a similar environment. The only physical (not numerical) 
approximation we make is linearization: our heat bath exchanges linear excitations (plane 
waves, or mesons) with the system, but not nonlinear ones (kinks, or baryons). Since we 
are only interested in pair creation and annihilation processes, the exchange of topological 
charge with the environment should have little impact on our results. It follows that we can 
approximately (the only approximation being linearization of the heat bath) identify our 
simulation time with the real physical one. The transition rates we measure have therefore 
direct physical meaning. 
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FIGURES 

FIG. 1. Cooling curves for a free-field system (m = 0.5, a = 1, L = 100) immersed in a simulated 
(squares) and real (solid line) zero-temperature heat baths. 

FIG. 2. Heating curves for a free-field system (m = 0.5, a = 1, L = 100) immersed in a simulated 
(squares) and real (solid line) heat baths. The equilibrium is reached at 9 = 1. 

FIG. 3. Autocorrelation functions of the endpoint (solid line) and midpoint (squares) fields for 
a free-field system (m = 0.5, a = 1, L = 100) at 9 = 1. 

FIG. 4. Temperature dependence of the kink number (logarithmic scale). The solid line is a fit 
for Meff = 0.737M. 

FIG. 5. Temperature dependence of the kink effective mass obtained by fitting pairs of consec- 
utive points from Figure ^ to Eq. ^ The inverse temperature /3 is that of a higher-temperature 
point in each pair. 



FIG. 6. The kink number autocorrelation time r at /? = 4 determined from Eq. 21, plotted as 
a function of the lag NAt. 

FIG. 7. Temperature dependence of the kink lifetime (logarithmic scale). The solid line is a fit 
for U = 1.06M. 
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FIGURE 5 
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